It would be great if @Andreas-Forster or somebody else from Shapemeans could have a look at it regarding the user interface and practical usability.

Thanks for taking a look:

• The diagonalize setting can be turned off such that the same (or very similar with withExtendedBasis = true) coefficients create the same shapes in the different models. The disadvantage is that the model then relies on a non-diagonalized basis, which some other provided functions don't expect. Calls to instance or sample are correct but marginalLikelihood provides the wrong result i think. i'll make that clearer in the documentation. I could also add an alternative function call that returns a matrix that could be used to move samples between the different spaces. Then i would see the argument to always diagonalize to avoid a problematic basis.
• to align only over translation withExtendedBasis should be set to false. Then there is also a one-to-one relationship between the model spaces. I'll make the documentation clearer.

I am not sure there is a suitable in-depth explanation of the used algorithm at the moment.

Thanks, @Andreas-Forster for testing it.
As we have correspondence, we could just recover the corresponding coefficients using the project method. To me it does not feel natural to return a projection matrix or keep the basis non-orthonormal.

What do you think?

I agree, keeping the basis non-orthonormal might confuse even more. And the projection matrix is not needed, as you mentioned, we have a project method.

I think updating the documentation and maybe later, if available link some more verbose document might be valuable.

I did some further testing about the influence of a non orthonormal basis in scalismo. It seems the DiscreteLowRankGaussian can handle such a non diagonalized basis. There are also various other operations already in the class that modify an orthonormal basis and lead to a non orthonormal shape model.

The only thing that i regularly do that shouldn't be done with a non orthonormal basis is directly accessing the variance field to check the overall variance of the model. without considering the basis this is incorrect for a non orthonormal model. All implementations i looked at are written to handle non orthonormal basis.

So i can see that this comes down to preference. Closer equivalence between model spaces versus a more intuitive formulation of the resulting model space.